計算機で一様擬似乱数を生成する.
実数 \(a < b\) に対し,区間 \([a, b]\) 上の 一様分布 (uniform distribution) とは,確率密度関数 \[p(x) = \begin{cases} 0 & \text{if} \ \ x < a \\ \frac{1}{b - a} & \text{if} \ \ a \leq x \leq b \\ 0 & \text{if} \ \ x > b \end{cases}\tag{2.1}\] を持つ確率分布.\(\mathcal{U}[a, b]\) と書く.区間を指定せずに一様分布というと,\(\mathcal{U}[0, 1]\) を指す.
set.seed(1) # seedの設定
runif(3) # 長さ3の一様擬似乱数
## [1] 0.2655087 0.3721239 0.5728534
関数 runif の出力は一様乱数と認識しても問題ないか?
\(\rightsquigarrow\) スペクトル (Spectral) 検定,コルモゴロフ・スミルノフ (Kolmogorov-Smirnov) 検定 を適用.
# スペクトル検定
set.seed(1)
x <- runif(1000)
y <- runif(1000)
z <- runif(1000)
plot_ly(x=x, y=y, z=z, type="scatter3d", mode="markers", size=0.1)
立方体に均等に分布しているように見える.
# コルモゴロフ・スミルノフ検定
set.seed(1)
u <- runif(1000)
ks.test(u, punif)
##
## One-sample Kolmogorov-Smirnov test
##
## data: u
## D = 0.024366, p-value = 0.5928
## alternative hypothesis: two-sided
一様分布から乖離しているのであれば,\(p\)-値は小さくなるはずであるが,ここでは一様擬似乱数を観測として見たときの \(p\)-値は小さくない.
\(\rightsquigarrow\) 一様擬似乱数 (初期設定はMersenne-Twister法) は一様分布からの独立な乱数列と見ても顕著な差が見られない.
\(a, b, y_0, n\): 事前に決める整数.
\[ y_m = (a y_{m - 1} + b) \mod n \ (m = 1, 2, \ldots) \\ x_m = \frac{y_m}{n}\] として\(\{ x_m: m = 1, 2, \ldots \}\) を出力する.
u <- numeric(3*1e3)
y <- 1234
n <- 2^31 - 1
a <- 65539
b <- 0
for(i in 1:length(u)){
y <- (a*y + b) %% n
u[i] <- y/n
}
u_mat <- matrix(u, nrow=3)
x <- u_mat[1,]
y <- u_mat[2,]
z <- u_mat[3,]
plot_ly(x=x, y=y, z=z, type="scatter3d", mode="markers", size=0.1)
2次元平面に整列している様子が見て取れる.
以下,理論的な記述の際は乱数を扱い,実際の計算では擬似乱数を用いるが,擬似乱数も乱数と表記する.
ここで,\(\mathcal{U}[0, 1]\) からの一様乱数が生成できるとする.\(X \sim \mathcal{U}[0, 1]\) に対して,線形変換
\[ Y = b + (a - b) X \] を施すことで,\(\mathcal{U}[a, b]\) に従う乱数が生成できる.
\(\because\) \(\mathcal{U}[a, b]\) は累積分布関数 \[ F(x) = \frac{x - a}{b - a} \ (a \leq x \leq b) \] を持つが,線形変換してできた \(Y\) の従う分布も \[ \mathbb{P}(Y \leq x) = \mathbb{P}\left( X \leq \frac{x - b}{a - b} \right) = \frac{x - b}{a - b} = F(x) \ (a \leq x \leq b) \] となり,同じ累積分布関数を持つため.
一様乱数をもとに,与えられた1次元の確率分布 \(P\) に従う乱数を生成する.
累積分布関数を \(F\) とすると, \[ \mathbb{P}(X \leq x) = F(x) \ (-\infty < x < \infty) \] となるような \(X\) を生成すればよい.
連続な確率分布であれば,任意の \(u \in (0, 1)\) に対して,\(F(x) = u\) となる \(x\) がただ1つに決まる.
\(\rightsquigarrow\) \(F^{-1}(u) = x\) で逆関数 \(F^{-1}\) が定義できる.このとき,一様乱数 \(U\) により, \[ X = F^{-1}(U) \] で \(P\) に従う確率変数を作れる.
\[\because \mathbb{P}(X \leq x) = \mathbb{P}(F^{-1}(U) \leq x) = \mathbb{P}(U \leq F(x)) = F(x) \] このようにして確率分布 \(P\) に従う乱数を生成する方法を 逆変換法 (inversion method) という.
\(\lambda > 0\) とする.指数分布 \(\mathcal{E}(\lambda)\) の累積分布関数は, \[ F(x) = 1 - e^{-\lambda x} \] である.したがって,一様乱数 \(U\) を用いて \[ F^{-1}(U) = - \frac{\log (1 - U)}{\lambda} \] によって指数乱数が発生できる.ここで,\(1 - U \sim \mathcal{U}[0, 1]\) より,\(-\log(U)/\lambda\) によっても指数乱数を生成できる.
R 言語で実装する.パラメータ \(\lambda\) を \(\lambda = 2.0\) とする.
set.seed(1)
lambda <- 2.0
u <- runif(3)
-log(u)/lambda
## [1] 0.6630539 0.4942642 0.2785628
組み込み関数 rexp でも生成する.
set.seed(1)
lambda <- 2.0
rexp(3, rate=lambda)
## [1] 0.37759092 0.59082139 0.07285336
同じ seed を用いても異なる値が出力されていることから,異なった実装がなされているようである.
計算時間の比較を行う.
lambda <- 2.0
n <- 1e5
microbenchmark(A <- -log(runif(n))/lambda, times=100)
## Unit: milliseconds
## expr min lq mean median uq
## A <- -log(runif(n))/lambda 2.932916 3.165062 3.313722 3.250917 3.386003
## max neval
## 5.024328 100
microbenchmark(B <- rexp(n, rate=lambda), times=100)
## Unit: milliseconds
## expr min lq mean median uq
## B <- rexp(n, rate = lambda) 3.465591 3.550562 3.755103 3.663531 3.852399
## max neval
## 5.572789 100
実行時間に大きな差はないようである.
実数 \(\mu\) と,正の実数 \(\sigma\) に対し,コーシー分布 (Cauchy distribution) \(\mathcal{C}(\mu, \sigma)\) は,確率密度関数 \[ p(x; \mu, \sigma) = \frac{1}{\pi \sigma \left( 1 + \frac{(x - \mu)^2}{\sigma^2} \right)} \ (-\infty < x < \infty) \] を持つ.\(\mu = 0, \sigma = 1\) としたコーシー分布 \(\mathcal{C}(0, 1)\) に注目する.累積分布関数は, \[ F(x) = \int_{-\infty}^x \frac{1}{\pi (1 + y^2)}dy \] となる.学部の解析学を思い出すと,この積分は, \[F(x) = \frac{\mathrm{arctan} x}{\pi} + \frac{1}{2}\] となる.よって, \[F^{-1}(u) = \tan \left( \pi u - \frac{\pi}{2} \right)\] を得る.したがって,一様乱数 \(U\) を用いて, \[X = \tan (\pi U - \pi/2)\] でコーシー乱数を生成できる.
また,コーシー分布は独立な2つの標準正規乱数 \(X_1, X_2\) の比 \[\frac{X_2}{X_1} \tag{2.2}\] の分布と等しい.以下で,逆変換法,(2.2)の方法,組み込み関数の3つの計算時間を比較する.
n <- 1e5
microbenchmark(A <- tan(pi*runif(n)-pi/2), times=100)
## Unit: milliseconds
## expr min lq mean median uq
## A <- tan(pi * runif(n) - pi/2) 4.955723 5.177404 5.373381 5.303579 5.505346
## max neval
## 6.806326 100
microbenchmark(B <- rnorm(n)/rnorm(n), times=100)
## Unit: milliseconds
## expr min lq mean median uq max
## B <- rnorm(n)/rnorm(n) 10.00244 10.32055 10.66322 10.46786 10.88738 12.83597
## neval
## 100
microbenchmark(C <- rcauchy(n), times=100)
## Unit: milliseconds
## expr min lq mean median uq max neval
## C <- rcauchy(n) 5.188717 5.419534 5.627272 5.552008 5.713074 7.21623 100
上述の逆変換法では,累積分布関数 \(F\) の逆関数 \(F^{-1}\) が存在することを仮定した.
離散の場合は逆関数を持たないが,一般に累積分布関数 \(F(x) = \mathbb{P}((-\infty, x])\) は次の性質を持つことを思い出す.
逆に,上記の性質を満たす関数 \(F\) は \(\mathbb{R}\) 上のある確率分布の累積分布関数である.
ここで,関数 \(F\) が単調非減少で右連続のとき, \[F^{-}(u) = \inf \{x \in \mathbb{R}: u \leq F(x) \}\] によって広義の逆関数 \(F^{-}: \mathbb{R} \to \mathbb{R}\) を定める.広義の逆関数 \(F^{-}\) に対しては,\(F^{-}(F(x)) = x\) や \(F(F^{-}(u)) = u\) は必ずしも成り立たないが,次が成り立つ.
\(x \in \mathbb{R}, u \in [0, 1]\) に対して, \[F^{-}(F(x)) \leq x\] \[ F(F^{-}(u)) \geq u \]
確率分布 \(P\) の累積分布関数 \(F\) の広義の逆関数 \(F^{-}\) を用いることで,連続の場合と同様に \(U \sim \mathcal{U}[0, 1]\) から \(F^{-}(U)\) を計算することで,\(P\) に従う乱数を生成できる.
有限集合 \(-\infty < x_1 < \cdots < x_n < \infty\) にのみ値をとる確率分布を考える.この分布の累積分布関数 \(F\) の広義の逆関数 \(F^{-}\) は, \[ F^{-}(u) = \begin{cases} x_1 & \text{if} \ \ 0 \leq u \leq F(x_1) \\ x_2 & \text{if} \ \ F(x_1) < u \leq F(x_2) \\ \vdots \\ x_n & \text{if} \ \ F(x_{n - 1}) < u \leq 1 \end{cases} \] である.可算無限集合でも同様.